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Abstract 



We present a novel numerical method to the time-harmonic inverse medium scattering 
problem of recovering the refractive index from near-field scattered data. The approach con- 
sists of two stages, one pruning step of detecting the scatterer support, and one resolution 
enhancing step with mixed regularization. The first step is strictly direct and of sampling 
type, and faithfully detects the scatterer support. The second step is an innovative appli- 
cation of nonsmooth mixed regularization, and it accurately resolves the scatterer sizes as 
well as intensities. The model is efficiently solved by a semi-smooth Newton-type method. 
Numerical results for two- and three-dimensional examples indicate that the approach is 
accurate, computationally efficient, and robust with respect to data noise. 
Key words: inverse medium scattering problem, reconstruction algorithm, sampling method, 
mixed regularization, semi-smooth Newton method 

1 Introduction 

In this work we study the inverse medium scattering problem (IMSP) of determining the re- 
fractive index from near-field measurements for time-harmonic wave propagation [ ]. Consider 
a homogeneous background space R d (d = 2, 3) that contains some inhomogeneous media oc- 
cupying a bounded domain Q. Let u inc — e lkx ' d be an incident plane wave, with the incident 
direction d E S d_1 and the wave number k. Then the total field u induced by the inhomogeneous 
medium scatterers satisfies the Helmholtz equation [8] 



where the function n{x) is the refractive index, i.e., the ratio of the wave speed in the homoge- 
neous background to that in the concerned medium at location x. The model describes not only 
time-harmonic acoustic wave propagation, but also electromagnetic wave propagation in either 
the transverse magnetic or transverse electric modes [ , Appendix]. 

Next we let 77 = (n 2 — l)/c 2 , which combines the relative refractive index n 2 — 1 with the wave 
number k. Clearly the coefficient 77 characterizes the material properties of the inhomogeneity 
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(i) 



and is supported in the scatterer Q C R d . We denote by / = rju the induced current, by G(x, y) 
the fundamental solution to the Helmholtz equation in the homogeneous background, i.e., 



G(x,y) 



l -H^(k\x-y\), d = 2, 

1 p ik\x-y\ 

d = 3, 



4tt \x — y\ ' 



where the function //"q refers to Hankel function of the first kind and zeroth-order. Then we can 
express the total field u as follows [ ] 

u = u™ c + [ G(x,y)I(y)dy. (2) 

By multiplying both sides of equation (2) by 77, we arrive at the following integral equation of 
the second kind for the induced current / 

I(x)= V u inc + V [ G(x,y)I(y)dy. (3) 

The reformulation (3) is numerically amenable since all computation is restricted to the scatterer 
support fi, which is much smaller than the whole space R d . Hence, the complexity is also very 
low. We will approximate the solution to (3) by the mid-point rule (cf. Appendix A). 

Next we let u s — u — u inc be the scattered field, which is measured on a closed curve/surface 
T enclosing the scatterers Then the IMSP is to retrieve the refractive index n 2 or equivalent ly 
the coefficient 77 from (possibly very noisy) measurements of the scattered field u s , corresponding 
to one or several incident fields. In the literature, a number of reconstruction techniques for the 
IMSP have been developed. These methods can be roughly divided into two groups: support 
detection and coefficient estimate. The former group (including MUSIC [ , ], linear sampling 
method [ , 3] and factorization method [ ] etc.) usually is of sampling type, and aims at detect- 
ing the scatterer support efficiently. The latter group generally relies on the idea of regularization 
(including Tikhonov regularization [2, 23], iterative regularization method [12, 1, 13], contrast 
source inversion [27], subspace regularization [ ] and propagation-backpropagation method [28]), 
and aims at retrieving a distributed estimate of the index function. These approaches generally 
are more expensive, but their results may profile the inhomogeneities more precisely. 

In this paper, we shall develop a novel two-stage numerical method for solving the IMSP. The 
first step employs a direct sampling method, recently developed in [16], to detect the scatterer 
support stably and accurately. It is based on the following index function 

( p) ||ui L2(r) ||G(,* p )|| L2(r) yXp(Eil > W 

where ft D Q is a sampling domain. Numerically, the method is strictly direct and does not incur 
any linear matrix operations. The method can detect reliably the scatterer support even in 
the presence of a large amount of data noises [ ]. In particular, a (much smaller) computational 
domain D C £1 can be determined from the index <E>, and furthermore the restriction $|^> may 
serve as a first approximation to the coefficient 77. 

The second step enhances the image resolution by a novel application of (nonsmooth) mixed 
regularization. With the approximation <E>|£> from the sampling step at hand, equation (3) gives 
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an approximate induced current / as well as an approximate total field u. Then we seek a 
regularized solution 77 to the linearized scattering equation 

/ G{x,y)u{y)r){y)dy = u s (x), (5) 

JD 

by an innovative regularization incorporating both L 1 and H 1 penalties. The L 1 penalty pro- 
motes the sparsity of the solution [25, 4, 11]. However, the estimate tends to be very spiky 
with the L 1 penalty used alone. Meanwhile, the conventional H 1 penalty can only yield globally 
smooth but often overly diffusive profiles. In this work we shall propose a novel mixed model that 
consists of a suitable combination of the L 1 and H 1 penalties. As we will see, this mixed model 
produces well clustered and yet distributed solutions, thereby overcoming the aforementioned 
drawbacks. It is the mixed model that enables us to obtain a clear and accurate reconstruction 
of the inclusions: The homogeneous background is vividly separated from the scatterers and 
both support and intensity of the inclusions are accurately resolved. 

Numerically, the L 1 penalty term gives rise to a nonsmooth optimality condition, which 
renders its direct numerical treatment inconvenient. Fortunately, by using complementarity 
functions, the optimality condition reduces to a coupled nonlinear system for the sought-for 
coefficient 77 and the Lagrangian multiplier, which is amenable to efficient numerical solution. 
We shall develop an efficient and stable semi-smooth Newton solver for the model via a primal- 
dual active-set strategy [ ]. Overall, the direct sampling method is very cheap and reduces 
the computational domain D in the mixed model (cf. (5)) significantly, which in turn makes 
the semi-smooth Newton method for the mixed model very fast. Hence, the proposed inverse 
scattering method is computationally very efficient. 

The rest of the paper is structured as follows. In Section 2, we recall a novel direct sampling 
method for screening the scatterer support fi, and derive thereby an initial guess to the coef- 
ficient 77. Then in Section 3 we develop an enhancement technique based on the idea of mixed 
regularization, and an efficient semi-smooth Newton solver. Finally, we present numerical re- 
sults for two- and three-dimensional examples to demonstrate the accuracy and efficiency of the 
proposed inverse scattering method. 



2 A direct sampling method 



In this section, we describe a direct sampling method to determine the shape of the scatterers, 
recently derived in [ ]. We only briefly recall the derivation, but refer the readers to [16] for 
more details. Consider a circular curve r (d — 2) or a spherical surface T (d = 3). Let G(x,x p ) 
be the fundamental solution in the homogeneous background. Then using the definitions of the 
fundamental solutions G(x, x p ) and G(x, x q ) and Green's second identity, we deduce 



2i$t(G(x p ,x q )) = 



G(x,x q 



^ dG(x, x p 
dn 



G{x,Xr 



dG(x,x q 
dn 



ds, 



(6) 



where the points x pi x q E f^r 5 the domain enclosed by the boundary T. 

Next we approximate the right hand side of identity (6) by means of the Sommerfeld radiation 
condition for the Helmholtz equation, i.e., 



dG(x 1 x p ) 
dn 



ikG(x, x p ) + h.o.t.. 
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Consequently, we arrive at the following approximate relation 

G(x,x p )G(x,x q )ds « k~ 1 ^s(G(x p , x q )) : 

which is valid if the points x p and Xq are not close to the boundary T. 

Now, we consider a sampling domain f2 enclosing the scatterer support fi. Upon dividing O 
into small elements {t^}, we may approximate the integral in the scattering relation (2) by 

u % x ) = L G(x,y)I(y)dy « y^WjGfayj), (7) 
Jn 

where the weight Wj is given by \rj\Ij with \rj\ being the volume of the jth element rj. The 
relation (7) is plausible if the induced current / is regular in each element and the elements {rj} 
are sufficiently fine. It also admits a nice physical interpretation: the scattered field u s at any 
fixed point x E T is a weighted average of that due to point scatterers located at {yj}- 
Combining the preceding two relations yields 

/ u s (x)G(x,x p )ds w Wj$s(G(yj, x p )). (8) 

Hence, if the sampling point x p is close to some point scatterer yj, i.e., yj E ri, then G(yj,x p ) 
is nearly singular and takes a very large value, contributing significantly to the sum in (8). 
Conversely, if the point x p is far away from all the physical point scatterers, then the sum will 
be very small due to the decay property of G(yj,x p ). 

These facts lead us to the index function <&(x p ) in (4) for any x p in the sampling region f2. 
In practice, if a point x p satisfies <&(x p ) ~ 1, then it likely lies within the support; whereas if 
$(xp) « 0, then the point x p most probably lies outside the support. Hence it serves as an 
indicator of the scatterer support. Consequently, we can determine a domain D C as one 
approximate scatterer support, and moreover, the restriction <f>\z> of the index $ to the domain 
D may be regarded as a first approximation to the sought-for coefficient r]. The subdomain D 
may be chosen as D = {x E Q : $(#) > /xmax^^ <E>(x)} with fi being a cut-off value, i.e., the 
union of elements whose index values are not less than a specified fraction of the largest index 
value over the sampling region £1. This determination of subdomain D will be adopted in our 
numerical experiments. 

This method is of sampling type (cf. [22] for an overview of existing sampling methods), and 
its flavor closely resembles multiple signal classification [24, 6, 10] and linear sampling method 
[7, 20]. However, unlike these existing techniques, it works with a few (e.g., one or two) incident 
waves, is highly tolerant with respect to noise, and involves only computing inner products with 
fundamental solutions rather than expensive matrix operations as in the other two techniques. 
The robustness of $ is attributed to the fact that the (high-frequency) noise is roughly orthogonal 
to the (smooth) fundamental solutions. 

3 Mixed regularization 

The direct sampling method in Section 2 extracts an accurate estimate D to the scatterer support 
as well as a reasonable initial guess to the medium coefficient 77, i.e., <&\d- In this part we 
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refine the approximation by exploiting the idea of nonsmooth mixed regularization. Given 
the approximation <&\d, we can compute the induced current / via (3) for each incident wave and 
the corresponding total field u on the domain D from (2). By substituting the approximation u 
into equation (2), we arrive at the following linearized problem 

y)u{y)r]{y)dy = u s (x) , x E Y. 

It is convenient to introduce a linear integral operator K : L 2 (D) \-> L 2 (F) defined by 

(Kr,)(x)= [ G(x,y)u(y)ri(y)dy. (9) 

JD 

We observe that the kernel G(x, y)u{y) is smooth due to the analyticity of fundamental 
functions G(x, y) away from the singularities and standard Sobolev smoothness of the total field 
u(y) (following from elliptic regularity theory [8]). Hence, the linear operator K : L 2 (D) \-> L 2 (T) 
is compact. As a consequence, the linearized problem (9) is ill-posed in the sense that small 
perturbations in the data can lead to huge deviations in the solution, according to the classical 
inverse theory [26], which is reminiscent of the severe ill-posedness of the IMSP, and its stable 
and accurate numerical solution calls for regularization techniques. 

We determine an enhanced estimate of the coefficient 77 from the linearized problem (9) by 
solving the following variational problem: 

min i f \Kr]-u s \ 2 ds + a [ \rj\dx + ^ [ \Vr]\ 2 dx. (10) 
2 Jv Jd 2 J D 

In comparison with more conventional regularization techniques, the most salient feature of the 
model (10) lies in two penalty terms: it contains both the L 1 penalty and the H 1 penalty, which 
exert drastically different a priori knowledge of the sought-for solution. The scalars a and /3 are 
regularization parameters controlling the strength of respective regularization. 

This variational problem (10) allows us to determine a coefficient 77 which is distributed 
yet clustered, i.e., exhibiting a clear groupwise sparsity structure in the canonical pixel basis. 
This a priori knowledge is plausible for localized inclusions/inhomogeneities in a homogeneous 
background. The model (10) is derived from the following widely accepted observations. The L 1 
penalty promotes the sparsity of the solution [25, 4, 11], i.e., the solution is very much localized. 
Hence, the estimated background is homogeneous. However, if the L 1 penalty is used alone, the 
solution tends to be very spiky and may miss numerous physically relevant pixels in the sought- 
for groups. That is, the desirable groupwise structure is not preserved. Meanwhile, the more 
conventional H 1 penalty [ ] yields a globally smooth profile, but the solution is often overly 
diffusive. Consequently, the overall structure stands out clearly, but the retrieved background 
is very blurry, which may lead to erroneous diagnosis of the number of the inclusions and their 
sizes. In order to preserve simultaneously these distinct features of the sought-for coefficient, i.e., 
sparsely distributed groupwise structures in a homogeneous background, a natural idea would 
be to combine the L 1 penalty with the H 1 penalty, in the hope of retaining the strengths of both 
models. As we shall see below, the idea does work very well, and the model is very effective for 
enhancing the resolution of the estimate to the coefficient 77. 

The general idea of mixed regularization, i.e., using multiple penalties, has proved very effec- 
tive in promoting several distinct features simultaneously. This general idea has been pursued 
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in the imaging community [17, 21]. However, to the best of our knowledge, the model (10) has 
not been explored in the literature, let alone its efficient and accurate numerical treatment. A 
detailed mathematical analysis of the model (10) is beyond the scope of the present paper. We 
refer interested readers to [ ] for some preliminary results on mixed regularization and to [ ] 
for a related model (elastic- net). 

To fully explore the potentials of the model (10), an efficient and accurate solver is required. 
We shall develop a semi-smooth Newton type method, which allows extracting very detailed 
features of the solutions to the model (10). The starting point of the algorithm is the necessary 
optimality condition of the variational problem (10), which reads 

K*Krj - f3Ar] - K*u 8 E -a^(r/), (11) 

where ij)(rj) — \\tjWli and the subdifferential dijj(r]) [18] is the set-valued signum function, which 
is defined pointwise as 

r i, if 7/0*0 >o, 

dtl>(v)(z) = { [-1,1], if^) = 0, 
[ -1, if 77(rc) < 0. 

Due to the convexity of the functional, the relation (11) is also a sufficient condition. Hence 
it suffices to solve the inclusion (11), for which there are several different ways, e.g., iterative soft 
shrinkage [9, 29], augmented Lagrangian method/alternating direction method or semi-smooth 
Newton method [18]. We shall develop a (new) semi-smooth Newton method to efficiently solve 
the inclusion (11). To this end, we first recall the complementarity condition [18] 

A = * + ,. (12) 
max(l, |A + cr]\) 

for any c > 0, which will be fixed at a constant in the implementation, and A serves as a Lagrange 
multiplier. It can be directly verified by pointwise inspection that the complementarity condition 
(12) is equivalent to the inclusion A E dtp(r]) (cf. [ ]). With the help of the complementarity 
condition (12), we arrive at the following equivalent nonlinear system in the primal variable 77 
and dual variable A: 

K*Kr] - /3Ar] - K*u 8 + a\ = 0, 

* ,it^ ii =°- 

max(l, |A + ct]\) 

Then we apply the semi-smooth Newton algorithm using a primal-dual active set strategy. 
The complete implementation is listed in Algorithm 1. The technical details for deriving the 
crucial Newton step (Step 5) are deferred to Appendix B, which involves damping and regular- 
ization. A natural choice of the stopping criterion at Step 6 is based on monitoring the change 
of the active set A = {x E D : |A + crj\ < 1}: if the active sets for two consecutive iterations 
coincide, then we can terminate the algorithm, cf. [18]. 

The main computational effort of the algorithm lies in the Newton update at Step 5: each 
iteration requires solving a (dense) linear system. We note that the dual variable A can be 
expressed in terms of the primal variable 77 and on the active set A, the coefficient 77 vanishes 
identically. Thus in practice, we solve only a linear system for 77 on the inactive set 1 = D \ A. 
An important feature of the algorithm is that the linear system becomes smaller and smaller 
and also less and less ill-conditioned as the iteration goes on, while the iterate captures more 
and more refined details of the nonhomogeneous medium regions. If the exact solution is indeed 
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Algorithm 1 Primal-dual active set method 
1: Initialize 77° and A , and set c > and k = 0. 
2: for k = 0, . . . , K do 

3: Set the active set A k and inactive set I k respectively by 

A k = {xeD: \X k + cr] k \ <1}, 
X k = {x e D : |A* + C77*| > 1}. 



4: Compute a and b by 



X k 1 7 \ k + CT] k 

and 6 



max(|A fe |,l) |A fe + c77 fe | 



and set d fc = |A fc + cr] k \ and F fc = afc*. 
5: Solve for (?7 fe+1 , A /c+1 ) from the system 



K*K V k+1 - $ A V k+1 - K*u s + a X k+1 = on X fc , 

0, 



rf fc -l v Jl max(|A fc |,l) 

rj k+1 = 0on/. 



6: Check the stopping criterion. 
7: end for 

8: output approximation rf K . 



sparse (many zero entries), then the system size, i.e., usually shrinks quickly as the iteration 
proceeds. The numerical experiments indicate that the convergence of the algorithm is rather 
steady and fast. 



4 Numerical experiments 

In this part, we present numerical results for several two- and three-dimensional examples to 
showcase the proposed two-stage inverse scattering method, for both exact and noisy data. The 
wave number k is fixed at 27T, and the wavelength is set to A = 1. The exact scattered field u s is 
obtained by first solving the integral equation (3) for the induced current / and then substituting 
the current / into the integral representation (2). Here the integral equation (3) is discretized 
by a mid-point rule; see Appendix A for details. The noisy scattered data u s s are generated 
pointwise by the formula 

u§(x) = u s {x) + e( max | u s (x) \ , 

where e refers to the relative noise level, and both the real and imaginary parts of the noise 
C — C( x ) follow the standard Gaussian distribution. The index <E>, its restriction and the 
enhanced approximation 77 by the mixed model will be displayed. As is mentioned in Section 2, 
we choose the subdomain D (approximate scatterer support) based on the formula D = {x E 
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(a) true scatterer (b) index <3> (c) index <!>\d (d) sparse recon. 

Figure 1: Numerical results for Example 1(a): (a) true scatterer, (b) index <3>, (c) index 
(restriction to the subdomain D) and (d) sparse reconstruction. The first and second rows refer 
to exact data and the data with 20% noise, respectively. 

/xmax ^ $(x)}, where the cut-off value /i is taken in the range (0.5, 0.7). The choice 
of the cutoff value /i affects directly the size of the domain D, but does not cause much effects 
on the reconstructions. 

Like in any regularization technique, an appropriate choice of regularization parameters 
(a, /3) in the mixed model (10) is crucial for the success of the proposed imaging algorithm. 
There have been a number of choice rules [15] for one single parameter, but very little is known 
about the mixed model. We shall choose the pair (a, f$) in a trial and error manner, which suffices 
our goal of illustrating the significant potentials of the mixed model for inverse scattering. In 
Algorithm 1, the parameter c is set to 50, and both rf and A are initialized to zero. The 
maximum number K of Newton iterations is 50. In all the experiments, the convergence of 
the algorithm is achieved within about 10 iterations. All the computations were performed on 
MATLAB 7.12.0 (R2011a) on a dual-core desktop computer with 2GB RAM. 

4.1 Two-dimensional examples 

Unless otherwise specified, one incident direction d is employed for two-dimensional problems, 
and it is fixed at 1) T - The scattered field u s is measured at 30 points uniformly distributed 

on a circle of radius 5. The sampling domain is fixed at [— 2, 2] 2 , which is divided into a uniform 
mesh consisting of small squares of width h — 0.01. The subdomain D for the integral equation 
(9) is divided into a coarser uniform mesh consisting of small squares of width 0.02. 
Our first example illustrates the method for two separate scatterers. 

Example 1. We consider two separate square scatterers in the following two scenarios 

(a) The scatterers are of width 0.2 and centered at (—0.8,-0.7) and (0.3,0.9), respectively, 
and the coefficient 77 in both region is 1. 
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Table 1: Regularization parameters (a, f3) for the examples. 



example 


1(a) 


1(b) 


2 


3 


e = 0% 
e = 20% 


(2.0e-6,1.5e-9) 
(3.0e-6,2.0e-9) 


(8.0e-6,1.4e-8) 
(8.5e-6,9.0e-9) 


(7.0e-6,1.0e-9) 
(7.0e-6,5.0e-9) 


(2.5e-9,4.0e-14) 
(2.5e-9,5.0e-14) 



(b) The scatterers are of width 0.3 and centered at (—0.25,0) and (0.25,0), respectively, and 
the coefficient 77 in the former and latter is 1.5 and 1, respectively. 

The two scatterers in Example 1(a) are well apart from each other. The recovery of the 
scatterer locations by the index $ is quite satisfactory, especially upon noting that we have 
just used one incident wave. Two distinct scatterers are observed for both exact data and the 
data with 20% noise, cf. Fig. 1(b). However, the magnitudes are inaccurate, and the estimate 
suffers from spurious oscillations in the homogeneous background, due to the ill-posed nature 
of the IMSP and the oscillating behavior of fundamental solutions. Nonetheless, two localized 
square subregions D (each of width 0.4) encompass the modes of the index <3>, see Fig. 1(c), and 
may be taken as an approximate scatterer support. In Fig. 1(c), the entire sampling domain 
Q is shown, and the two small squares represent the approximate support D. Outside of the 
domain D, the index is set to zero, i.e., identical with homogeneous background, and will 
not be updated during the enhancing step. The enhancing step is initialized with and the 
results are shown in Fig. 1(d). The regularization parameters for getting the reconstructions, 
which are determined in a trial-and-error manner, are presented in Table 1. The enhancement 
of the approximation over the domain D is significant: the recovered background is now 
mostly homogeneous, and the magnitudes and sizes of the recovered scatterers agree well with 
the exact ones. This shows clearly the significant potentials of the proposed mixed regularization 
for inverse scattering problems. The numbers in Table 1 also sheds valuable insights into the 
mixed model (10): the value of the regularization parameter a is much larger than that of /3. 
Hence, the L 1 penalty plays a predominant role in ensuring the sparsity of the solution, whereas 
the H 1 penalty yields a locally smooth structure. 

The two scatterers in Example 1(b) stay very close to each other, and thus it is rather 
challenging for precise numerical reconstruction. The detection of the scatterer locations by 
the index <E>, see Fig. 2(b), is still very impressive. In particular, it clearly distinguishes two 
separate scatterers with their locations correctly retrieved, and this remains stable for data with 
up to 20% noise. The mixed regularization is performed on a square D (of width 1) enclosing 
the two modes in the index see Fig. 2(c). This choice of the inversion domain D allows 
possibly connecting of the modes. However, the enhancement correctly recognizes two separate 
scatterers, with their magnitudes and sizes in excellent agreement with the exact ones. Also the 
estimated background is very crispy. Surprisingly, the estimate deteriorates only slightly in that 
the right scatterer elongates a little bit towards the left scatterer as the noise level e increases 
from to 20%. Although not presented, we would like to note that for this particular example, 
the reconstructions are still reasonable for data with 30% noise. Hence the proposed inverse 
scattering method is very robust with respect to the data noise. 

Next we consider a ring-shaped scatterer. 

Example 2. The scatterer is one ring-shaped square located at the origin, with the outer and 
inner side lengths being 0.6 and OA, respectively. The coefficient 77 of the scatterer is 1. Two 
incident directions d\ — -4=(1, 1) T and d^ — -4=(1, — 1) T are considered. 
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(a) true scatterer 



(b) index <E> 



(c) index 



(d) sparse recon. 



Figure 2: Numerical results for Example 1(b): (a) true scatterer, (b) index <3>, (c) index <fr\r> 
(restriction to the subdomain D) and (d) sparse reconstruction. From top to bottom, the rows 
refer to exact data, noisy data with 10% and 20%, respectively. 

Ring-shaped scatterer represents one of most challenging objects to recover, and it is highly 
nontrivial even with multiple scattered field data sets, especially noting the ring has a small 
thickness. It has been observed that one single incident field is insufficient to completely resolve 
the ring structure, and only some parts of the ring can be resolved, depending on the incident 
direction d [ ]. Hence we employ two incident waves in the directions d\ — ^(1,1) T and 

^2 = ^(1,— 1) T in order to yield sufficient amount of information about the scatterer, and 
accordingly, the index function <E> is defined as follows 

i 

where the function <3>^ refers to the index for the ith data set. The numerical results with the 
exact data and 20% noise in the data are shown in Fig. 3. With just two incident waves, 
the index <3> can provide a quite reasonable estimate of the ring shape. Despite some small 
oscillations, the overall profile stands out clearly, and remains very stable for up to 20% noise 
in the data. The enhancing step via mixed regularization provides a very crispy estimate of the 
ring structure: the recovered scatterer has a clear ring structure, which agrees excellently with 
the exact one in terms of both magnitude and size. The presence of 20% data noise causes visible 
deterioration to the reconstruction, see Fig. 3(d). Nonetheless, the enhanced reconstruction still 
exhibits a clear ring shape, and it represents a very good approximation to the true scatterer 
upon noting the large amount of data noise. 
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(a) true scatterer (b) index <3> (c) index <!>\d (d) sparse recon. 

Figure 3: Numerical results for Example 2: (a) true scatterer, (b) index <E>, (c) index 
(restriction to the subdomain D) and (d) sparse reconstruction. The first and second rows refer 
to exact data and the data with 20% noise, respectively. 

4.2 Three-dimensional example 

Our last example shows the feasibility of the method for three-dimensional problems. 

Example 3. We consider two cubic scatterers of width 0.1 centered at (0.35,0.15,0.15) and 
(—0.35, 0.15, 0.15), respectively. One single incident field with direction d— ^(1, 1, 1) T is used, 
and the coefficient r\ of the scatterers is taken to be 1. 

The scattered field u s is measured at 600 points uniformly distributed on the surface Y of a 
cubic of width 5, (i.e., 10 points in each direction). To simulate the scattered field data, we take 
the sampling domain ft to be the cubic [— 1, l] 3 , which is divided into a uniform mesh consisting 
of small cubes of width h = 0.01. The inversion domain D for the integral equation (9) is divided 
into a coarser mesh consisting of small cubes of width 0.03. 

The numerical results for Example 3 with exact data are shown in Fig. 4(b), where each row 
represents a cross-sectional image along the second coordinate axis X2- The scatterer support 
estimated by the index $ agrees reasonably with the exact one, and away from the boundary 
of the true scatterers, the magnitude of $ decreases quickly. However, the reconstructed profile 
is slightly diffusive in comparison with the exact one, which is reminiscent of the decay prop- 
erty of fundamental solutions. The nonsmooth mixed regularization (10) is carried out on two 
cubic subregions (of width 0.36A), cf. Fig. 4(c). Like before, a significant improvement in the 
resolution is observed: the sparse estimate is much more localized in comparison with the index 
<3>, and also the magnitude is close to the exact one; see Fig. 4(d). The presence of 20% data 
noise does not worsen much the index <3> and the sparse reconstruction, cf. Fig. 5. Hence the 
reconstruction algorithm is highly tolerant with respect to data noise. 

Lastly, we briefly comment on the computational efficiency of the overall procedure. The 
first step with the index involves only computing inner products and is embarrassingly cheap 
and easily parallelized. The accuracy of the support detection is quite satisfactory, and thus a 
large portion of the sampling domain ft can be pruned from inversion, i.e., \D\ <C \ft\. Hence, 
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the enhancement via mixed regularization is also rather efficient. 

5 Concluding remarks 

We have presented a novel two-stage inverse scattering method for the inverse medium scattering 
problem of recovering the refractive index from near-field scattered data. The efficiency and 
accuracy of the method stem from accurate support detection by the sampling strategy and 
group sparsity-promoting of the mixed regularization technique. The former is computationally 
very efficient, and reduces greatly the computational domain for the more expensive inversion 
via nonsmooth mixed regularization, while the latter achieves an enhanced resolution with the 
magnitudes and sizes comparable with the exact ones. The numerical results for two- and 
three-dimensional examples clearly confirm these observations. 

These promising experimental results raise a number of interesting questions for further 
studies. First, the potentials of mixed regularization have been clearly demonstrated. It is of 
great interest to shed theoretical insights into the model as well as to design efficient acceler- 
ation strategies, which for three-dimensional problems remains very challenging. Some partial 
theoretical results can be found in [14]. Also of much practical relevance is an automated choice 
of regularization parameters. Second, the reconstructions were obtained with the linearized 
model, which represents only an approximation to the genuine nonlinear IMSP model. It would 
be interesting to justify the excellent performance of the linearization procedure. Third, the ro- 
bustness of the approach to noise is outstanding when compared with more conventional inverse 
scattering algorithms, especially noting the limited data for inversion. The mechanism of the 
robustness is not yet clear. 
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A Numerical method for forward scattering 

We denote by J the index set of grid points of a uniformly distributed mesh with a mesh size 
h > and consider square cells 

Bj — Bjij2 ~ ( x ji' x j 2 ) + I - \ ' f] x [f ' |] 

for every tuple j = (ji, J2) belonging to the index set J. Assume that the domain U jG j Bj contains 
the scatterer support fi. We use the mid-point quadrature rule to evaluate the operator K, and 
hence the integral (3) is approximated by 

h-VkYl Gkjtjh 2 = Vk u inc (x k ) 
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(a) true scatterer (b) index <3> (c) index $|£> (d) sparse recon. 

Figure 4: Numerical results for Example 3 with exact data: (a) true scatter, (b) index <E>, (c) 
index <f>\p (restriction to the subdomain D) and (d) sparse reconstruction. From the top to 
bottom: the cross sectional images at X2 — 1.07, 1.10, 1.13, 1.16, 1.19, and 1.22, respectively. 



13 




(a) true scatterer (b) index <3> (c) index $|£> (d) sparse recon. 

Figure 5: Numerical results for Example 3 with 20% noise in the data: (a) true scatter, (b) index 
<3>, (c) index <f>\]j (restriction to the subdomain D) and (d) sparse reconstruction. From the top 
to bottom: the cross sectional images at X2 — 1.07, 1.10, 1.13, 1.16, 1.19, and 1.22 respectively. 
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where Ik = I{xk) an d % = v( x k), an d the off-diagonal entries Gkj and the diagonal entries 
are given by Gkj = G(xk,xj) and 




respectively. The diagonal entries can be accurately computed by tensor-product Gaussian 
quadrature rules. The resulting system can be solved using standard numerical solvers, e.g., 
Gaussian elimination, if the cardinality of the index set J is medium, and iterative solvers like 
GMRES. The extension of the procedure to 3D problems is straightforward. 



B Semi-smooth Newton method 

In this part, we derive a semi-smooth Newton method for minimizing (10). The optimality 
condition of the variational problem reads 

( K*Kr] + aX- (3Ar] - K*u 8 = 0, 
max(l, |A + crj\) 

where A is the Lagrange multiplier (dual variable). The second line, the complementarity func- 
tion, equivalently expresses the inclusion A E ^H^Hl 1 ? which can be checked directly by pointwise 
inspection. Thereby, we effectively transforms the inclusion (11) into a numerically amenable 
nonlinear system. It follows directly from the complementarity relation 

A = * + (13) 
max(l, |A + C7]\) 

that on the active set A = {x E D : |A + crj\(x) < 1}, 77 vanishes identically. Otherwise, both 
the dual variable A and the primal variable 77 need to be solved. We shall solve the system by a 
semi-smooth Newton method [ ] . First observe that the Newton step (with the increments for 
A and 77 denoted by 5 A and £77, respectively) applied to the following reformulation of equation 
(13) (on the set 1 = D \ A) 

A|A + C77I — A + C77 = 

is given by 

I A + C7]\5X + X .^ + Cr, . \SX + cStj] - (5X + cSr]) + AIA + cq\ - (A + cry) = 0, 
IA + C77I 

or equivalently with the notation A + = A + 5X and 77+ = 77 + £77, we have 

A + |A + cr]\ + A [A + + cr] + ] = A|A + crj\ + [A + + cr] + }. 
Next we apply the idea of damping and regularization to the equation and thus get 
A + |A + cq\ + 6[X + + cy ? + ] ,^ + Cyy , ^— — = [A + + ct? + ] + 9\X + cq\ X 



\X + C77I max(|A|, 1) max(|A|, 1) 
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Here, the purpose of the regularization step max ^| -n is to automatically constrain the dual 
variable A to [—1,1]. The damping factor 9 is automatically selected to achieve the stability. To 
this end, we let d — |A + crj], rj = d — 1, a = max (j A | fj ? and 6 = p^p^jy- We arrive at 

\ + (rj + 1) + (9[A + + cr/ + ]a6 = [A + + cr? + ] + 6>ad. 

Thus we have 

1 f)d 

A+ = ^r-pr-ril - 9ab}c V + + ^-^a 
r] + 6ab l 11 r] + 6ab 

To arrive at a simple iteration scheme, we set ~+Q ab = 1, i.e., = < 1. Consequently, we 
obtain a simple iteration 

■ 1 — ab , A 

d — 1 max(|A|, 1) 

where we have used the relation = ^E^- Substituting this into the first equation gives 

K*Kr] + + ac^-^V - /3Ar/+ = K*u s - a (14) 

d — 1 max(|A|, 1) 

We note that one only needs to solve equation (14) on the inactive set X, since on the active 
set A, there always holds r] + = 0. This has an enormous computational consequence: the size 
of the linear system in (14) can be very small if \I\ is small, i.e., the solution is sparse. This 
last relation shows also clearly the sparsity of the solution, and this provides a crispy estimate 
of the background. Upon obtaining the solution 77+ , one can update A + on the sets I and A 
according to the second and the first equation, respectively. Lastly, we would like to remark on 
the consistency of the scheme: if the sequence generated by the semi-smooth Newton method 
converges, then the limit satisfies the complementarity relation (13) as desired. 
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